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ABSTRACT 

We construct phase-space distribution functions for the oblate, cuspy mass models 
of Sridhar & Touma, which may contain a central point mass (black hole) and have 
potentials of Stackel form in parabolic coordinates. The density in the ST models is 
proportional to a power r~'^ of the radius, with < 7 < 1. We derive distribution 
functions f{E,Lz) for the scale-free ST models (no black hole) using a power series 
of the energy E and the component Lz of the angular momentum parallel to the 
symmetry axis. We use the contour integral method of Hunter & Qian to construct 
f{E,Lz) for ST models with central black holes, and employ the scheme introduced 
by Dejonghe & de Zeeuw to derive more general distribution functions which depend 
on E, Lz and the exact third integral I3. We find that self-consistent two- and three- 
integral distribution functions exist for all values < 7 < 1. 

Key words: celestial mechanics - stellar dynamics - galaxies: kinematics and dy- 
namics - galaxies: structure - galaxies: central black holes. 



1 INTRODUCTION 

Observations of the nuclei of elliptical galaxies using the 
Hubble Space Telescope have revealed that the luminos- 
ity density diverges towards the centre as a power r ' of 



Lauer et al. 1995; CaroUo et 



the radius (Jaffe et al. 1994 

|al. 1997t [Rest et al. 2001] ). Man y and perhaps all of the 
nuclei host a central black hole ( Richstone et al. 199!: ; |d^ 
|Zeeuw 2001 ). Dynamical models for such cusped galaxies 
include the axisymmetric power-law galaxies introduced by 
Toomre (1982) and Evans (1994), which have spheroidal 
potentials and simple f{E,Lz) distribution functions, with 
E the orbital energy and the component of the angu- 
lar momentum parallel to the symmetry axis. Similar mod- 
els with spheroidal densities were constructed by Dehnen 
& Gerhard (1994) and Qian et al. (1995, hereafter Q95). 
Most orbits in these models admit an approximate third in- 
tegral of motion. Self-consistent three-integral distribution 
functions were constructed for some special scale-free cases 
by Schwarzschild's (1979) numerical orbit superposition 
method (Richstone 1980, 1984; Levison & Richstone 1985) 
and by analytic means (de Zeeuw, Evans & Schwarzschild 
1996; Evans, Hafner & de Zeeuw 1997, hereafter EHZ). 

Q95 constructed two-integral distribution functions 
f{E,Lz) for axisymmetric power-law spheroidal densities 
containing a central black hole by means of the contour 
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integral method of Hunter & Qian (1993, hereafter HQ). 
They showed that self-consistent f{E,Lz)'s exist for oblate 
scale-free spheroids with < 7 < 3, but found an upper 
bound for the axial ratios of self-consistent prolate models 
of this kind. Their study also revealed that the presence of 
a central black hole limits the region of parameter space for 
which physical distribution functions (i.e., f{E,Lz) > 0) 
exist to 7 > 1/2 (for all axis ratios). De Bruijne et al. 
(1996) studied oblate spheroidal cusps in the radial range 
where the spherical potential of the black hole dominates, 
and constructed analytic constant-anisotropy three-integral 
distribution functions f{E,Lz,L), with L the modulus of 
the angular momentum. These distribution functions remain 
physical also when < 7 < 1/2, and suggest that shallow- 
cusped axisymmetric galaxies with central black holes can 
in fact be constructed, a result confirmed by application of 
Schwarzschild's numerical method (Verolme, priv. comm.). 
Self-consistent axisymmetric models with central black holes 
designed to fit the observed photometry and kinematics 
of specific galaxies were constructed with Schwarzschild's 
method by, e.g., van der Marel et al. (1998), Cretton et al. 
(1999, 2000) and Gebhardt et al. (2000). 

The oblate models introduced by Sridhar & Touma 
(1997a, hereafter ST) provide the only set of cuspy mod- 
els with a central black hole for which all orbits have an 
exact third integral of motion Is, so that, in principle, exact 
distribution functions can be constructed. The potential of 
these models is of Stackel form in parabolic coordinates, and 
this causes the density to be significantly fiattened, with a 
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shape that is fixed once the cusp slope is chosen, and a total 
mass that is infinite. The oblate axisymmetric models with 
a potential of S tackel form in p rolate spheroidal coordinates 
(Kuzmin 1956; ie Zeeuw 1985) also admit an exact third in- 
tegral and include a large set of models with a range of den- 
sity profiles and shapes with finite total mass. However, they 
all have constant density cores and can not contain a central 



point mass without destroying the separability (de Zeeuw 



Peletiei & Franx 1986). Oblate density distributions with 



potentials that are separable in spherical coordinates have 
densities that become negative at large distances (Lynden- 
Bell 1962b) . Here we study self-consistent distribution func- 
tions for the oblate ST models. We show that, by contrast to 
the cusped spheroidal densities of Q95, the ST models with 
a central black hole do have consistent f{E, Lz) distribution 
functions for all values of the cusp slope < 7 < 1. 

In we summarize the properties of the ST models. 
Without a black hole, the ST models are scale-free. We dis- 
cuss their distribution functions in §^ and investigate the 
general case in We summarize our conclusions in §0. 



2 SRIDHAR-TOUMA MODELS 

A comprehensive description of the mass models and the 
orbit structure can be found in ST. Here we collect the rele- 
vant properties, derive the fundamental integral equation for 
the distribution function, and close with a brief discussion 
of the velocity moments. 



2.1 Potential, density and orbits 

The motion in the axisymmetric ST models separates in 
parabolic coordinates rj) in the meridional plane. We de- 
fine them as 



5 — r(l + cos ( 



r(l-cos0), C,'7>0, 



(1) 



where r is the polar radius and 6 is the co-latitude. Our 
definition of rj differs from that in ST by an overall sign, 
which removes the need to use I?;! in many expressions. In 
these coordinates the gravitational potential of an ST model 
can be written as 



2if + n'^-^') - 2GM 



K>0, 



(2) 



where K is a positive constant, G is the gravitational con- 
stant, and M denotes the mass of a central black hole (point 
mass). The density distribution p{£,,rf) associated with the 
potential (|) follows from eq. (18) of ST. When M = 0, 
the potential and the density are proportional to r'^'"' and 
r respectively, i.e., they are scale-free, but the density is 
non- negative only for < 7 < 1. The equipotential surfaces 
are approximately spheroidal. Their axis ratio b/a (defined 
by the condition V{Q,h) = V{a,Q)) equals 1/2 for all values 
of 7. As a result, the surfaces of constant density are dim- 
pled along the short (2) axis, and the dimple deepens, i.e., 
the density distribution becomes increasingly toroidal, with 
increasing 7. 

For subsequent use, we record here the expressions for 
the ST potential-density pairs in terms of the standard 



spherical polar coordinates (r, 9, 

GM 
r ' 

p{r,e) = poT~^S{e) + M5{r), 



V{r,e) = Kr^P{e) 



(3) 



where po ~ 2(3 — 7)(1 — 'y)K/nG, S is the Dirac delta- 
function, and 



P{e) = (l + cos( 



\3—l 



+ (1 — cos ( 



\3-~f 



S(e) = (2-7-cos6')(l+cos6i)^"^ 

-h(2-7 + cos6l)(l-cos6')^"^. 

In terms of cylindrical polar coordinates (cc7, tf} 

p{^,z) = g{[(2-7)r-2](r + z)^-^ 

+ [(2 - 7)r + z](r-z)'-^}+M(5(r), 



(4) 



(5) 



where r = V tu^ + 2^. In what follows, we take units such 
that K = \ and po = 1- 

The separation of the Hamilton- Ja cobi equat ion results 
in a third integral of motion given by dPars 196E| ) 

+ ^ - - GM, (6) 



73 = 2ip\ + 2e 
or, equivalently 



3-7 



/3 = -27?p^ - 2rf-'' 



7-2 

^ + Er] + GM, 



(7) 



with and pr, the momenta conjugate to ^ and ry, respec- 
tively. Here E denotes the orbital energy and = p^ = 
r^4> sin^ 9 the component of the angular momentum parallel 
to the symmetry axis, both of which are integrals of mo- 
tion as well. All orbits with 7^ are short-axis tubes, 
which may be either symmetric or asymmetric with respect 
to the equatorial plane. Those with = are confined to a 
plane of constant (j), and are centrophilic when a black hole is 
present. The orbital families are illustrated in Figure 2 of ST. 



2.2 Distribution functions 

In integrable systems, the distribution function depends on 
the phase-space coordinates through the isolating integrals 
of motion (Jeans 1915; Lynden-Bell 1962a), which for ax- 
isymmetric systems means / = f{E, L^, I3). The mass den- 
sity p is related to / through the integral 



p{r,9) = f{E,L,,h)A' 



8 



f{E,L,,h 



d{^,V,<t))d{Pi,Pr,,P4,) 



d{r,9,4>) d{E,h,L, 



dEdLAh, (8) 



where v = {vr,V0,v^)'^ is the velocity vector in the spheri- 
cal coordinates, and the factor of eight occurs because both 
E and I3 are quadratic in the velocities, and we consider 
only distribution functions that are even in L^, so that 
f{E,—Lz,I-i) = /(EjLzjIs). The Jacobians in the integral 
(g) follow from expressions (Q), ^ and and are given 
by 



d{r,9,^) 



2r sin 9, 
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d{E,h,L,) 







Asme^[h-I^][I+-h] 



(9) 



where 



lt{E,L,,ri) 



3-7 



2ri 



+ Eri + GM 



We define the equivalent surface density E(r, 9) by 
E(r,e') = 2rsin6lp(r,6i), 
and substitute expression in the integral to obtain 



(10) 

(11) 



di. 



d/s' 



V{r,e) 

where 



(12) 



(13) 



This is the fundamental integral equation for the distribu- 
tion functions of the ST models. According to ( |l^ ) , the dis- 
tribution functions obtained from ( |l^ will be even in L^. 
We assign only positive values to I/zmaxi which means that 
all stellar orbits have a definite sense of rotation around the 
axis of symmetry. 



-(2C-fr,){«|)]-p|^, 

-{2n + i){vl)]~p^. (18) 

The relations with the familiar second moments in spherical 
coordinates are 



pi^e) = [np{4) + ^p{-"v)] . 

p{vrVo) = ^ - P{4)] ■ (19) 

These relations are valid for any oblate density p in a po- 
tential V that is separable in parabolic coordinates. In the 
two-integral limit, we have (u^) = {vg) and {vrVo) — 0, or, 
equivalently, (y^) = {v^). 

Using (Fm) one can show that the second velocity mo- 
ments for the self-consistent ST models are all infinite. The 
problem occurs at S = oo. For the case Ai = it also follows 
through the application of eq. (2.14) of EHZ, which gives a 
divergent integral for the ^-dependence of the stresses, and 
confirms that this property is shared by all weak cusps with 
< 7 < 1. Scale-free steep cusps with 7 > 1 do not have 
this problem, because their density falls ofi^ sufficiently fast 
with radius. 



2.3 Velocity moments 

The velocity moments p{v\^v^v^) are defined by 



/ I m n\ III I m n r / j- \ i.3 



(14) 



where the integrations are carried out over all possible veloc- 
ities. The velocity components are related to the momenta 
in parabolic coordinates through 

p^ = C'i = Cv(, pr, = M^rj = Mvn, p,j,=J^^<j> = J^Vij,, (15) 
where C, M and N are the metric coefficients defined as 

Transforming to the integration variables E, Lz and I3 gives 
(cf. eq. Q) 

p{4<^;;)=2^+^e'^^' 
f{E,L,,h)L^ 



■ dhdL.dE, 



. . ^ .... (17) 

(13-/3-)— (/3+-/3) — 

where the integration limits are identical to those of the 
integral (|l|). 

We consider only the second moments. It is not diffi- 
cult to show that the velocity ellipsoid is aligned with the 
parabolic coordinate system, so that (w^Wrj) = {v^v^} = 
{"VriV^} = 0. This leaves three non-zero elements of the stress 
tensor, {v^), and (w^), which are connected to the po- 

tential V and the density p by two Jeans equations, 

dpivj) ^ p 



3 SCALE-FREE ST MODELS 

In the absence of a central black hole, the ST models are 
scale-free, and it is natural to consider distribution func- 
tions which are also scale-free. This restriction simplifies the 
fundamental integral equation, as we show in this section. 



3.1 Fundamental integral equation 

When Ad — 0, the integrals and I3 can be scaled by 
the energy E. We consider scale-free distribution functions 
f(E,Lz,h) of the form 



f{E,L,,l3)=E^g{E''l3,E'L,), 



(20) 



where the exponents p, q and s are real numbers. We define 
the dimensionless variables 

u=^, v^EiE", w = L,E\ (21) 
from which we obtain 

dE = -E—, dl3 = E-''dv, dL,=E~'dw. (22) 
u 

We choose the values of p, q and s so that the radius r factors 
out of the integral equation (^2|). This occurs when; 
10 1 11 

(23) 



2 2-7 2-7- 2-7 2 

It follows that scale-free distributions for the ST models are 
of the form (po|), with p, q and s given in (p3|). 

In terms of the new variables (^), equation (^^ be- 
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(a) (6) (c) 

Figure 1. (a) The variation of R{9) = cr{0) — crM{d) for 7 = 0.7 and N = 28. (b) The distribution function Q{w) for < wJ = ui/uimax < 1- 
(c) The history of the objective function versus N. J has converged to its minimum at 10~^^. The numerical accuracy is saturated 
because of truncation and round-off errors. 



(j(ff) — u '^-1 Au dw dv 



(24) 



-V- V+-V 



v^{u,w,0) 

where the integration limits are given by 

sin 9 



'W+ 



(u,e) = ^2(i-it)u5^ 



v+(u, w, 



(1 — cos6')u2-T 



P(6l)-2w(l-cos( 



x2-7 



2 ^ p{e)^ 

-W U 2-T 



V-{u,W,9) — 



{1 + cos 9)u-^-~< 



3 — r 

P{9)'^ 



2(l-cose»)2. 
2tt(l+cos6l)^"^-P(6l) 



2 2_ pre)^ 



2(l+cose')2 



and the left hand side is defined as 
25(61) sine 



o{9)^ 



P{9)^ 



(25) 



(26) 



Hence, the problem of finding f{E, L^, I3) is reduced to solv- 
ing equation (M) for Q{v,w). 



3.2 Fricke series for flEjL^) 

We first consider the (even part of the) two-integral distri- 
bution function /{EjL^), which is determined uniquely by 
the density distribution. One way to obtain it is to assume 
that Q only depends on w. In this case equation (Eih reads 



1 iu_|_(u,9) 

a{9) — n I u 2-7 du I C/(w)dw. 







(27) 



This integral equation must be solved for G{w) > 0. We 
assume the solution in the form of power series (Fricke 1952) 



(28) 



where we have to determine an {n — 0, ... , N) and we choose 
A*' based on the required accuracy of the solutions. We sub- 
stitute the series ( psj ) into the integral equation ^t\) and 
obtain 



JV 

O-Jv(0) = angn{9), 
n=0 



(29) 



where the functions gn{9) are given by (Gradshteyn & 
Ryzhik 1980) 

1 w^(u,e) 

f ^7~3 f n 

gn{9) — n I u 2-T du / w^dw 



B,J±2,!1±2 , (30) 



TT / V2stti9 



n + l 



P{9Y 



and B is the Beta-function. We determine the a„ from the 
requirement that on{9) converges to (j{9) as TV is increased. 
We think of mean convergence and therefore attempt to min- 
imize the objective function 



/ \ ['j{9)-aN{e)fd9, 



(31) 



with respect to the variations of the coefficients a„. This 
is the well-known method of Bubnov-Galerkin (e.g., Reddy 
1986). The function J has a local extremum if 

-i^ = 0, j = 0,...,iV. (32) 

Therefore, we are left with a set of linear algebraic equations 
for a„ as 



K a = b. 



(33) 



where the matrix K = [kij] is the so-called stiffness matrix, 
a — (ao, ai, . . . , ajv)^ is the vector of unknowns and b — 
(60, bi, • • • , biv)^ is a constant vector. The elements of K 
and b are given by 



hj = k,, = J g^{9)gj{e)d9, b, = J g,{9)ai9)d9. (34) 
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We note that 0(6) = cr(7r — 9) because of the symmetry with 
respect to the equatorial plane. Since the Qnid) have this 
property as well, n can take both even and odd values. The 
convergence to o-{6) is guaranteed if limiv^oo ^ = 0. 

As an example, we set 7 = 0.7 and follow the procedure 
mentioned above. We compute a„ for different values of A'' 
and increase N until J attains its minimum (« 10"^'^) and 
the accuracy of the numerical solutions is saturated. Our 
convergence condition is satisfied for N — 28. We find that 
aN{0) agrees with aiO) to an accuracy of lO"'^, which is con- 
sistent with the minimum value of J' . The series expansion 
of Q{'w) is also convergent in the domain < u) < Wmax 
where we have defined 



w+{- 



(35) 



so that lUmax correspouds to the circular orbit in the equa- 
torial plane. Figure la shows the variation of the residual 
function R{6) — cr{9) — (jm(Q) versus 9. The corresponding 
distribution function, which is a well-defined positive func- 
tion, has been plotted versus w = w/Wina.:^ and shown in 
Figure 16. By studying the history of J versus A'^ (Figure 
Ic) it follows that the extremum of J is indeed a mini- 
mum. Furthermore, the envelope of the residual function is 
almost uniform, which indicates optimal fitting of ct{0). Our 
numerical experiments show that the speed of convergence 
increases when 7 — > 0. 

An alternative method for the determination of the a„ 
coefficients would be to expand a {9) and gn{9) in Fourier 
series, and to calculate a„ by comparing the coefficients of 
the sine functions on both sides of (p^). However, the eval- 
uation of the Fourier coefficients considerably increases the 
required computational effort, and we have not followed this 
approach. Direct evaluation of f{E, L^) by means of the HQ 
method is discussed in 

Our results show that, just as for the scale-free 
spheroids of Q95 and the scale-free power-law galaxies of 
Evans (1994), the scale- free ST models admit self-consistent 
f(E,Lz) distribution functions. In the former models the 
function G{w) is monotonic, but for the ST models it has a 
maximum. 



3.3 Distribution functions f{E,Lz,Is) 

Dejonghe & de Zeeuw (1988, hereafter DZ) constructed 
three-integral distribution functions for Kuzmin's (1956) 
model using a generalized Fricke's (1952) method. They ex- 
panded the given model density p in terms of the p oten- 
tial V and the polar cylindrical radius zaj = ^yx^ + y^, and 
wrote the distribution function / as the sum of two parts: 
/ — /2-I-/3. /2 and /a are two- and three- integral distribu- 
tion functions, respectively. One can integrate fs to obtain 
the corresponding density pa, which is subtracted from p. 
The remaining density p — pa is then reproduced by /2. DZ 
computed the coefficients of the Fricke expansion by direct 
comparison of the series representations for p and the inte- 
gral of /. 

We can construct three-integral distribution functions 
by solving (|2^) for Q{v,w) by means of a method devel- 
oped by DZ for axisymmetric systems with Stackel poten- 
tials. This can be considered as a perturbative approach in 



which one exploits the existence of two-integral distribution 
functions. We assume Q has the form 



Q{v,w) = Go{v,w) +gi{w). 



(36) 



Then, we choose a specific form for Goi'Vj'w) and compute 
the resulting density a{9). Subtracting 5-{9) from (y{9) leaves 
a residual function R{9). Finally, we determine the two- 
integral Gi{w) that is consistent with R{9). 

We first consider simple monomial forms for Qo, 

Qo{v,w) ^ ao,mnv"^w". (37) 

Any combinations of ( ^ ) can also be chosen. By substi- 
tuting ( ^tI ) in the fundamental integral equation (|24[), and 
carrying out the integrations over v, w and u, we find 



a{9) = ao^rnngO,mn{9). 

The explicit expression for go,-, 



(38) 



i{0) is (see Appendix A) 



m m — 1 m — i — j i 



go,; 



'W-EE E E(t)(70M(0 



1=0 j=0 fe=0 i=0 



x(-l) 



m — i — j — k-\-l i-)i + fc + - 



^B(i + i i) 



xB (^+k+^, + 
{l + n + 2j+2l)P{9)^ 

where 

P,(e) = i±^, P,i9) ''''^ 



P{9)^ ' 



Piey- 



(39) 



(40) 



We now calculate R{9) = a{9) — a{9) and solve the following 
equation for Gi{'w) 



R{9) 



Qi{w)dw. 



(41) 



The solution of this equation for Qi{w) can be inserted in 
( p6| ) to determine Q{v,w). The performance of this technique 
depends on the initial choice of Qoi^jW). 

As an example, we construct a three-integral distribu- 
tion function for the same case discussed in §3.2. We take 
7 = 0.7, ao,22 = 50 00 and assume m — n = 2. This gives us 
the basis function 50,22 (S) that is plotted in Figure 2a. The 
residual function of this step, R{9), is displayed in Figure 
2b. We now assume Qi{w) = X^^q b^w^ and find K so that 
R{6) is reproduced with the best available accuracy, i.e.. 



{u,6) 



R{e) ^ R{e) 



27-3 
u ^--1 du 



w'^dw. 



(42) 



The constants are calculated by minimizing the objective 
function Ji = ^[R{e) - R(9)]'^d9 with respect to the 
variations of bk. Numerical computation shows that con- 
verges to a minimum of 5.3 x lO"^'^ by taking K = 15. The 
global residual function. Raid) — R{9) — R{9), is plotted in 
Figure 2 c. The three-integral distribution function obtained 
from ( |36| ) is positive for all possible values of u and w. The 
residual functions of the two- and three-integral distribu- 
tion functions, displayed in Figures la and 2c, have simi- 
lar behaviours. Their envelopes are nearly identical, and the 
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2E-06 



0.3 0.6 0.9 1.2 1.5 




(a) 



Figure 2. (a) The basis function go,22{d) for 7 = 0.7 and m = n = 2. (b) The residual function R{8) when cr{6) is approximated by 
Qo{v,w). (c) The global residual function Rq{9). 



maximum deviation from cr{9) is approximately equal to the 
square root of the minimum value of the objective function, 
as expected. The minimum value of the objective function 
is the limit of available numerical accuracy. 



4 ST MODELS WITH CENTRAL BLACK 
HOLES 

In the presence of a central black hole, i.e., when M 7^ in 
equation (^), the ST models still have separable potentials 
but lose their scale-freeness. The construction of distribution 
functions by means of the series solutions of §^ becomes more 
complicated, even for the two-integral case. We overcome 
this problem by using the contour integral method of Hunter 
& Qian (1993). 

4.1 Construction of /(iJ, L^) 

The expressions for the density and potential of the ST 
models in cylindrical polar coordinates (ro, <jf>, z) are given 
in equation 1^. We consider the relative potential = —V 
corresponding to the relative energy £ — — ^v^ (Binney & 
Tremaine 1987) with v the modulus of the velocity vector. 

The potential 'i' can be split into 'I' = ^'bh + ^* where 
^BH ~ GM/r and \['* is the potential induced by the density 
(^). The value of a point (zo, z) is obtained through 

Ll = 2uj' [*(^^ z')-E] = -2v. f'^^ff\ ^=^^, (43) 

where Wc is the radius of the circular orbit of energy £ in 
the equatorial plane, which is obtained from 



dro2 



(44) 



At a fixed energy £, Lz takes its maximum value for zu = -ujc 
and 2 = 0. We denote this maximum by Lc = Lzizoc) and 
define w = L^/Lc, which is equal to ui/iUmax of equation 
(H). For the ST models equation (^ reads 



(4 - 7)^^-^ + £7^7, 



iGM : 



0. 



(45) 



It is convenient to express (|45D in terms of the dimension- 
less parameter ( = [*BH(ti7f7o)/|**(^, 0)|] > 0. Thus, we 
obtain GM = 2C'cu'^~''' and equation ( [45| ) is equivalent to 



-£/(4-7-C)]~. 



(46) 



It follows that C > 1 corresponds to the black hole sphere 
of infiuence, i.e., the region where the potential of the black 
hole dominates, and < < 1 outside this region. 

To apply the HQ method, we follow Q95 and express the 
density p in terms of ^ and zo'^ as p{'^,-co'^) by eliminating 
between the expressions for the density and potential 
given in equation (fel). The contour integral solution of HQ 
for axisymmetric models has the form 



47r2iv^ 



.{£)+] 



2(p-£) 



dp 



(47) 



which is evaluated in the complex p-plane. Following HQ 
and Q95, we use the subscript 1 to indicate a partial deriva- 
tive with respect to the first argument of p. The symbol 
[*I'env(f )+] indicates the contour of integration, which starts 
from on the lower side of p-plane, crosses the real axis 
at and ends at on the upper side of p-plane. 

For the ST models we have '^00 = —00 and ^'env(f ) = 
'^{zol, 0) = 2ro^~^(C - 1)- The value of *cnv belongs to the 
interval [^'min, ^max] of physically achievable potentials on 
the real axis. To determine /5ii(p,ti7^) we use the implicit 
relation 



Pii(p, ^ 



/ 2 2\ 
P22{j^ , Z ) 



P2(n7^, Z^)*22(tIJ^ 
[^2{vJ^,z'^)f 



(48) 



where each subscript 2 indicates partial differentiation with 
respect to the second argument, z^ . For a given pair 
{p,vj^ — L\l\2(p — f)]} on the contour of integration, the 
variable z? is supplied to (^) through solving the following 
nonlinear complex equation 



2(p-£)' 



P = 0. 



(49) 



We use the modified Newton method for solving (|49|), and 
use the recursive formula 



2 

2n+l 



V'{z'^ 



di' 



for obtaining the (n -I- l)th estimate of the root of T'{z^). 
The coefficient 1/2-' allows for a univariate search along the 
gradient of V with the aim of taking the best step size to- 
wards the final answer. The integer exponent j is determined 
through 
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www 
(a) (6) (c) 

Figure 3. Physical distribution functions constructed for ST models using the HQ contour integral method. In all cases we have set 
S = —I. Panel a with C = corresponds to scale-free models in the absence of central black holes. Panel b corresponds to = 0.3 outside 
the black hole sphere of influence. Panel c with = 1.5 corresponds to the regions inside the black hole sphere of influence. 



j = mm{fc ; < fc < fcn,ax, ||P(4 + )|| < 117^(4)11}, (51) 

which guarantees a uniform and stable convergence. This 
algorithm fails if fc > fcmax and we have to change the initial 
choice for starting the recursion (^). 

An appropriate contour is the one used by Q95, and 
defined by 



P=*cnv(f)+/ 



1 — sec 



(I) 



+i/isini9, -TT < 1? < TT, (52) 



with / and h positive constants. The maximum width of 
the integration contour is controlled by h, while I adjusts 
the location of the local maximum of the contour. The sin- 
gularities (poles and branch points) of pii(p, tu'^) play an 
important role in the integration along [\['env(f )+]• Unfor- 
tunately, due to the implicit evaluation of pii(p,r*7^), we 
have no clear idea about the possible singularities except 
for the branch point p = £ that corresponds to -^p — £ in 
the denominator of the integrand of (kTI). To gain a better 
sense, we changed the width of our contour and investigated 
the existence of singularities by monitoring the value of dis- 
tribution function. We set I = h = c|^'env(f )| and evaluated 
the integral ^ for c = 0.2, 0.5, 1, 5 and 10. The results 
were the same in all cases indicating that p — £ is the only 
singular point on the real axis and the integrand does not 
have any complex conjugate singularities. Our computations 
show that larger values of c give more accurate results, and 
therefore we have used c = 5 throughout. 

We evaluate the contour integral (^^ using a Gaussian 
quadrature. We carry out a change of independent variable 
as A = cosi9 (0 < I? < tt) and obtain 



-I ^ +ih\ 



\/2(l + A) 
which transforms (^) to 



-dA 



f{£,L 
where 



2< (a+i&) — (a — ib) b 



-1 < A < +1, (53) 



(54) 



a + ib = / 5(A) 



dA 



ff(A) 



-L—zz hlAlA 

^(1 + A) 



Pii pW, 



Li 



2[p(A)-£-], 



(55) 



Physical distribution functions correspond to & > 0. Assum- 
ing W{X) = l/vT^-A^ as the weight function, one can apply 
the A'^-point Gauss-Chebyshev quadrature formula (Press et 
al. 1992) and obtain 



5(A) 



dA 



where 



— A 

A' 



7r(j-l/2) 

N 



(56) 



(57) 



We solve (|5C| ) with an accuracy of 10" and increase 
A until an accuracy of 10~* is obtained in the integration 
of dss] ) using ([s^. For instance, we constructed distribution 
functions for £ = —1 and several values of C, and 7. Figure 
3 shows the results for C, = Q, 0.3 and 1.5. Scale-free models 
without central black holes correspond to ^ = 0. As can be 
seen in Figure 3a, the graph of 7 = 0.7 is in agreement with 
Q{w) of Figure 16, which was constructed by means of the 
Fricke series. 

By increasing the value of C, (moving into the black hole 
sphere of influence), the maxima of f{E, Lz)'s are shifted to 
the high-Lz orbits (Figure 3). This means that sufficiently 
close to the central black hole, more nearly circular orbits 
are needed to maintain self-consistency in the ST models. 



4.2 Distribution functions f{E,Lz,Is) 

Three-integral distribution functions of ST models with cen- 
tral black holes can be derived as in §^ for the scale- free mod- 
els. We again write /(£;, Lz, /a) = f2{E, L^) + ef3{E, L^, I3), 
assume simple forms for /a and determine the density pa — 
e J fgd^v. We then subtract p3 from p and generate /2 by 
the residual density P2 ~ p ~ P3 using the method of HQ. 

For example, we assume /s = E~'^I^5{Lz) where 5 is the 
Dirac delta function. The corresponding density becomes 

for models with central 



(this is a consequence of eq. 



resp 

m 
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Figure 4. The two-integral /2 part of tlie three-integral distribution functions of ST models for E = —2. The left and right panels 
correspond to the scale-free ST models and ST models with central black holes, respectively. The value of f = 0.5 is relevant to the 
regions outside the black hole sphere of influence. 



black holes) 

P3 



2^ 



\/3 



V 



1/2 



where 

ci = -2tj'-^ + GM, C2 = 2^ 



GM, 



(58) 
(59) 



with ^ and rj the parabolic coordinates. The potential func- 
tion V is given by (H). It is now straightforward to construct 
f2{S,Lz) iiom p2{uj = /9(t47^, 2^) — p3(ci7^, 2^) by follow- 
ing the procedure in §4.1 . We have set e — 0.02, £ — —2 and 
have constructed /2(£,iz) > for several choices of 7 and 
C,. The results have been demonstrated in Figure 4. As this 
figure shows, f2{S, Lz) is oscillatory and steeply increases as 
«J — > 0. Our three-integral distribution functions exist both 
in the absence and in the presence of a central black hole 
for all values of < 7 < 1. 



5 DISCUSSION 

We have constructed two- and three-integral distribution 
functions for the cusped oblate models with central black 
holes introduced by Sridhar & Touma (1997a). We have 
computed f{E,Lz) for the case without a black hole by 
means of a Fricke (1952) series, and have confirmed the re- 
sult by means of the contour integral method of Hunter & 
Qian (1993). The distribution function f{E,Lz) is of the 
form E^Qiw), with p — — i — ^'^'^ ^ — Lz/Lc where 
Lc{E) is the value of Lz for the circular orbit of energy E. 
These distribution functions are non-negative for all values 
of the cusp slope < 7 < 1, in agreement with the results 
obtained by Evans (1994) for scale-free power-law galaxies, 
and by Q95 for scale-free oblate spheroidal densities. By 
contrast to these other models, the function G{w) is not 
monotonic for the strongly dimpled shape of the ST models. 

The ST models with a central black hole are not scale- 
free, but the HQ contour method shows that in this case 
again self-consistent f{E, Lz) exist. This makes the ST mod- 
els significantly different from spheroidal cusps with central 



black holes, which do not admit a self-consistent f{E,Lz) 
when 7 < 1/2 (Q95; de Bruijne et al. 1996). We ascribe 
this difference to the dimpled shape of the ST models, as 
the density distribution can be considered as the weighted 
integral of two-integral components S{E — E*)5{Lz — LI), 
each of which have toroidal shapes. This result is unlikely to 
depend on the separability of the ST models, or on the de- 
tails of the orbit structure, as the computation of f{E,Lz) 
does not require any knowledge of the orbits, or indeed of 
the existence of an exact third integral. We speculate that 
the power- law galaxies of Evans (1994) can have a physi- 
cal f{E, Lz) distribution function even when a central black 
hole is added. As the density p of these models is a simple 
function of and vj^ (see Appendix D of Evans & de Zeeuw 
1994), these distribution functions can be found by m eans 



of the HQ method following the procedure described in §4.1 

Nevertheless, the simple and exact form of I3 in the ST 
models makes it possible to construct distribution functions 
f[E,Lz,H) for these models, by means of a scheme intro- 
duced by Dejonghe & de Zeeuw (1988). These three-integral 
ST models have stars on the symmetric short-axis tube or- 
bits as well as on pairs of reflected banana orbits (to ensure 
the symmetry with respect to the equatorial plane). The ST 
models have special shapes, but have a somewhat richer dy- 
namical structure than, e.g., the oblate models with Stackel 
potentials in spheroidal coordinates, which contain only one 
major orbit family, the short-axis tubes. 

The two-dimensional versions of the ST models (Sridhar 
& Touma 1997b) are elongated discs, and no self-consistent 
distribution functions exist (Syer & Zhao 1998). The non- 
self-consistency of these models is related to the nature of 
the banana orbits that deposit much mass far from the major 
axis where the density is maximum. Although the meridional 
motions of the oblate ST models suffer this shortcoming too, 
their measure in the [E, Lz, I3) space is zero, and the short- 
axis tubes with Lz ^ have a sufficient variety of shapes to 
allow a range of self-consistent distribution functions. 

Scale-free models with shallow cusps have mass distri- 
butions that diverge strongly at large radii, and infinite pro- 
jected surface densities and stress tensors. Realistic models 
clearly require a radial profile that falls off more steeply at 
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larger radii. Q95 and de Bruijne et al. (1996) have shown 
that at small and large radii the scale-free models provide 
insight into the dynamics of these more general models. The 
ST models similarly provide insight into the nature of galac- 
tic nuclei with a central black hole and a shallow luminosity 
cusp. In particular they indicate that the detailed shape of 
the surfaces of constant density may have a significant in- 
fluence on the nature of the stellar velocity distribution. 
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APPENDIX A: DETERMINATION OF BASIS 
FUNCTIONS 

The functions go.mn(d) are given by 



gO,mniO)= Ju ^-~< du J dw J 



v_ {u,w,6) 



(Al) 



where V- and v+ are given in equation (|2^). Define the new 
variable 



< ^ < 1, 



(A2) 



and rewrite ( [A1| ) as 



go,^4o)-iu-^du L'^d«;/i!^^±i^;±=^=M::ii^. (as) 



v^^Ci-m) 



The inner integral in (A3) can be evaluated by writing out 



the binomial expansion for [v- + (w+ — W-)/i]™, and using 
the definition of the beta-function. The result is 



m 

^(7)^;— («+-«_)'B(^^-i,i; 



(A4) 



We now collect the w-dependent terms of V- and v+ — V- 
and write 

V- =vi{u, 9)+w%2{u, 9), V+—V- =V3{u, 9)+w%4{u, 6),{A5) 



where 
vi{u, 9) 



1 + cos 9 r„ „x2 , 

^[2u(l+cos6»)" ^- 



P{9)~ 



-P{9)]u~ , 



P{9)~ 



V2{U,9) = U 2-7 

2(l + cos6') 
V3{u,9) = 2u^ {l-u)P{9)~^ , 



V4{u,9) — 



j_P(6l) — 



(A6) 



By combining these results, and carrying out the integration 
over w, we obtain 



/IIL 1 1 L — L L 
.^d.5:5:5:(7)(7')(;) 
^ n — (\ A — A? — n 



1=0 i=o ;=o 



xB(j -I- i,i)«r~'~'' 



n + l + 2j + 2l 



n + l + 2j + 21 



(A7) 



© 0000 RAS, MNRAS 000, 000-000 



10 M. A. Jalali and P. T. de Zeeuw 

Finally, we replace u™"'^-' with its binomial expansion, col- 
lect the ii-dependent terms and then integrate over u. This 
gives us equation (691). 
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